/*******************************************************************************
	- Simulation of performance changes by the end of century (robustness checks)
*******************************************************************************/

********************************************************************************
* initial settings
********************************************************************************

version 14.2
clear 
drop _all
set more off

*----------------------------------

// set your own local path where you put the JAERE_Submission folder
if c(username) == "zhenx" {
	global LOCALPATH "C:/Users/zhenx/Dropbox/Research/Temperature and Athletes/temp_acclimatization/Code" 
}
global PATH "$LOCALPATH/JAERE_Submission"

********************************************************************************
* Table A6: estimated performance change under different temperature scenarios
********************************************************************************

global RESULT "$PATH/Results/Simulations/Robustness"
cd "$PATH/Results/Simulations/Baseline"

foreach scenario in "RCP45" "1F" "2F" "4F" {

	*** load data set 
		use "$PATH/Data/TF_Data.dta", clear

	*** keep only endurance events
		keep if event_type == "Endurance"

	*** drop training temperature lower than 40F to mitigate noise
		/*avoid the case where athletes might get trained indoors during extreme cold weathers*/
		drop if Team_tavg <= 40 

	*** set scenario
		if "`scenario'" == "RCP45" {
			gen Meet_GDDPtemp = Meet_RCP45temp
			gen Team_GDDPtemp = Team_RCP45temp 
		}
		if "`scenario'" == "1F" {
			gen Meet_GDDPtemp = Meet_tavg + 1
			gen Team_GDDPtemp = Team_tavg + 1
		}
		if "`scenario'" == "2F" {
			gen Meet_GDDPtemp = Meet_tavg + 2
			gen Team_GDDPtemp = Team_tavg + 2
		}
		if "`scenario'" == "4F" {
			gen Meet_GDDPtemp = Meet_tavg + 4
			gen Team_GDDPtemp = Team_tavg + 4
		}

	*** only keep weather variables
		keep athlete Meet_tavg Team_tavg Meet_GDDPtemp Team_GDDPtemp Team_zip norm_result

	*** set other control variables at the mean level
		gen wind = -0.0000445
		gen Meet_preciptotal = 0.0819
	    gen Meet_ozone = 0.0367131
	    gen travel_dist = 277.2341

	*** drop observations with missing data
		egen rowmiss = rowmiss(Meet_tavg Team_tavg Meet_GDDPtemp Team_GDDPtemp)
		drop if rowmiss > 0

	*** construct temperature variables
		* team temperature 
		mkspline team_c_linspl1 50 team_c_linspl2 60 team_c_linspl3 = Team_tavg, di
		mkspline team_f_linspl1 50 team_f_linspl2 60 team_f_linspl3 = Team_GDDPtemp, di

		* meet temperature - cubic spline
		mkspline2 meet_c_cubspl= Meet_tavg, cubic di knots(49.6 56.7 63.8)
		mkspline2 meet_f_cubspl= Meet_GDDPtemp, cubic di knots(49.6 56.7 63.8)

		* meet temperature - linear spline with 2 knots
		mkspline meet_c_linspl1 50 meet_c_linspl2 60 meet_c_linspl3 = Meet_tavg, di
		mkspline meet_f_linspl1 50 meet_f_linspl2 60 meet_f_linspl3 = Meet_GDDPtemp, di

		* meet temperature - linear spline with 3 knots
		mkspline meet_c_3linspl1 50 meet_c_3linspl2 60 meet_c_3linspl3 70 meet_c_3linspl4 = Meet_tavg, di
		mkspline meet_f_3linspl1 50 meet_f_3linspl2 60 meet_f_3linspl3 70 meet_f_3linspl4 = Meet_GDDPtemp, di

		* backup meet and team temperature
		gen Meet_tavg_bkp = Meet_tavg 
		gen Team_tavg_bkp = Team_tavg 

	*** initialize variables to store linear splines of team pre-exposure temperatures
		forvalues i = 1/3 {
			gen team_linspl`i' = .
		}

	*** linear spline with 3 knots
		forvalues i = 1/4 {
			gen meet_3linspl`i' = .
		}

		* scenario 1
		estimates use "Semi_Parametric_Linear_Spline_3knots_Homo"

		forvalues i = 1/4 {
			replace meet_3linspl`i' = meet_c_3linspl`i'
		}
		predict yhat_c1_linear_spline3, xb

		forvalues i = 1/4 {
			replace meet_3linspl`i' = meet_f_3linspl`i'
		}
		predict yhat_f1_linear_spline3, xb

		* scenario 2
		estimates use "Semi_Parametric_Linear_Spline_3knots_Hetero"

		forvalues i = 1/4 {
			replace meet_3linspl`i' = meet_c_3linspl`i'
		}
		forvalues i = 1/3 {
			replace team_linspl`i' = team_c_linspl`i'
		}
		predict yhat_c2_linear_spline3, xb

		forvalues i = 1/4 {
			replace meet_3linspl`i' = meet_f_3linspl`i'
		}
		forvalues i = 1/3 {
			replace team_linspl`i' = team_f_linspl`i'
		}
		predict yhat_f2_linear_spline3, xb


	*** quadratic
		* scenario 1
		estimates use "Semi_Parametric_Quadratic_Homo"

		replace Meet_tavg = Meet_tavg_bkp
		predict yhat_c1_quadratic, xb

		replace Meet_tavg = Meet_GDDPtemp
		predict yhat_f1_quadratic, xb

		* scenario 2
		estimates use "Semi_Parametric_Quadratic_Hetero"
		replace Meet_tavg = Meet_tavg_bkp
		forvalues i = 1/3 {
			replace team_linspl`i' = team_c_linspl`i'
		}
		predict yhat_c2_quadratic, xb

		replace Meet_tavg = Meet_GDDPtemp
		forvalues i = 1/3 {
			replace team_linspl`i' = team_f_linspl`i'
		}
		predict yhat_f2_quadratic, xb

	*** cubic
		* scenario 1
		estimates use "Semi_Parametric_Cubic_Homo"

		replace Meet_tavg = Meet_tavg_bkp
		predict yhat_c1_cubic, xb

		replace Meet_tavg = Meet_GDDPtemp
		predict yhat_f1_cubic, xb

		* scenario 2
		estimates use "Semi_Parametric_Cubic_Hetero"
		replace Meet_tavg = Meet_tavg_bkp
		forvalues i = 1/3 {
			replace team_linspl`i' = team_c_linspl`i'
		}
		predict yhat_c2_cubic, xb

		replace Meet_tavg = Meet_GDDPtemp
		forvalues i = 1/3 {
			replace team_linspl`i' = team_f_linspl`i'
		}
		predict yhat_f2_cubic, xb

	*** linear spline with 2 knots
		forvalues i = 1/3 {
			gen meet_linspl`i' = .
		}

		* scenario 1
		estimates use "Semi_Parametric_Linear_Spline_2Knots_Homo"
		forvalues i = 1/3 {
			replace meet_linspl`i' = meet_c_linspl`i'
		}
		predict yhat_c1_linear_spline2, xb

		forvalues i = 1/3 {
			replace meet_linspl`i' = meet_f_linspl`i'
		}
		predict yhat_f1_linear_spline2, xb

		* scenario 2
		estimates use "Semi_Parametric_Linear_Spline_2Knots_Hetero"
		forvalues i = 1/3 {
			replace meet_linspl`i' = meet_c_linspl`i'
		}
		forvalues i = 1/3 {
			replace team_linspl`i' = team_c_linspl`i'
		}
		predict yhat_c2_linear_spline2, xb

		forvalues i = 1/3 {
			replace meet_linspl`i' = meet_f_linspl`i'
		}
		forvalues i = 1/3 {
			replace team_linspl`i' = team_f_linspl`i'
		}
		predict yhat_f2_linear_spline2, xb

	*** cubic spline
		forvalues i = 1/2 {
			gen meet_cubspl`i' = .
		}

		* scenario 1
		estimates use "Semi_Parametric_Cubic_Spline_Homo"
		forvalues i = 1/2 {
			replace meet_cubspl`i' = meet_c_cubspl`i'
		}
		predict yhat_c1_cubic_spline, xb

		forvalues i = 1/2 {
			replace meet_cubspl`i' = meet_f_cubspl`i'
		}
		predict yhat_f1_cubic_spline, xb

		* scenario 2
		estimates use "Semi_Parametric_Cubic_Spline_Hetero"
		forvalues i = 1/2 {
			replace meet_cubspl`i' = meet_c_cubspl`i'
		}
		forvalues i = 1/3 {
			replace team_linspl`i' = team_c_linspl`i'
		}
		predict yhat_c2_cubic_spline, xb

		forvalues i = 1/2 {
			replace meet_cubspl`i' = meet_f_cubspl`i'
		}
		forvalues i = 1/3 {
			replace team_linspl`i' = team_f_linspl`i'
		}
		predict yhat_f2_cubic_spline, xb

	*** calculate percentage changes
		foreach model in linear_spline3 quadratic cubic linear_spline2 cubic_spline {
			gen percent1_`model' = (yhat_f1_`model' - yhat_c1_`model') / abs(norm_result)
			gen percent2_`model' = (yhat_f2_`model' - yhat_c2_`model') / abs(norm_result)
		}

	*** calculate mean changes
		collapse (mean) percent1* percent2*, by(ath Team_zip Team_tavg)
		foreach model in linear_spline3 quadratic cubic linear_spline2 cubic_spline {
		forvalues i = 1/2 {
			replace percent`i'_`model' = percent`i'_`model' * 100
		}
		}

		collapse (mean)	percent1_linear_spline3 percent1_cubic_spline percent1_linear_spline2 percent1_quadratic percent1_cubic ///
						percent2_linear_spline3 percent2_cubic_spline percent2_linear_spline2 percent2_quadratic percent2_cubic ///
				 (semean) se1percent_linear_spline3 = percent1_linear_spline3 ///
				 		  se1percent_cubic_spline   = percent1_cubic_spline ///
				 		  se1percent_linear_spline2 = percent1_linear_spline2 ///
				 		  se1percent_quadratic 	 	= percent1_quadratic ///
				 		  se1percent_cubic  		= percent1_cubic ///
				 		  se2percent_linear_spline3 = percent2_linear_spline3 ///
				 		  se2percent_cubic_spline   = percent2_cubic_spline ///
				 		  se2percent_linear_spline2 = percent2_linear_spline2 ///
				 		  se2percent_quadratic      = percent2_quadratic ///
				 		  se2percent_cubic          = percent2_cubic 
		gen i = _n
		reshape long percent1 percent2 se1percent se2percent, i(i) j(model) string
		order percent1 se1percent percent2 se2percent
		replace i = 2 if model == "_quadratic"
		replace i = 3 if model == "_cubic"
		replace i = 4 if model == "_linear_spline2"
		replace i = 5 if model == "_cubic_spline"
		sort i

	*** export results 
		export excel using "$RESULT/Estimated_Performance_Change_`scenario'.xlsx", firstrow(var) replace

}


********************************************************************************
* Table A7: estimated performance change using average of minimum and maximum of
*           historical temperature
********************************************************************************

*** set results folder
	cd "$PATH/Results/Simulations/Robustness/Min Max Avg"

*----------------------------regression estimates------------------------------*

*** load data set 
	use "$PATH/Data/TF_Data.dta", clear

*** define control variables and fixed effects
	global weather_control wind Meet_preciptotal Meet_ozone travel_dist 
	global absorb_var athlete season_week event_code##(venue home1) year_code i.Meet_dew_bin

*** set clustered standard error
	global cluster_se meet_code athlete

*** drop the observations that have missing values on key variables
	drop if Meet_tavg == . | Team_tavg == . | Meet_preciptotal ==. | Meet_ozone == .

*** drop training temperature lower than 40F to mitigate noise
	/*avoid the case where athletes might get trained indoors during extreme cold weathers*/
	drop if Team_tavg <= 40 

*** use average of tmin and tmax
	replace Meet_tavg = (Meet_tmax + Meet_tmin) / 2
	replace Team_tavg = (Team_tmax + Team_tmin_7) / 2

*** keep only endurance events
	keep if event_type == "Endurance"

*** construct linear splines for pre-exposure temperatures
	mkspline team_linspl1 50 team_linspl2 60 team_linspl3 = Team_tavg, di

*** linear spline with 3 knots
	cap drop meet_3linspl*
	mkspline meet_3linspl1 50 meet_3linspl2 60 meet_3linspl3 70 meet_3linspl4  = Meet_tavg, di

	reghdfe norm_result c.meet_3linspl* $weather_control, absorb($absorb_var) vce(cluster meet_code athlete)
	estimates save "Semi_Parametric_Linear_Spline_3knots_Homo", replace

	reghdfe norm_result c.meet_3linspl* (c.team_linspl1 c.team_linspl2 c.team_linspl3)#c.meet_3linspl* $weather_control, absorb($absorb_var) vce(cluster meet_code athlete)
	estimates save "Semi_Parametric_Linear_Spline_3knots_Hetero", replace

*** quadratic
	reghdfe norm_result Meet_tavg c.Meet_tavg#c.Meet_tavg $weather_control, absorb($absorb_var) vce(cluster meet_code athlete)
	estimates save "Semi_Parametric_Quadratic_Homo", replace

	reghdfe norm_result Meet_tavg c.Meet_tavg#c.Meet_tavg (c.team_linspl1 c.team_linspl2 c.team_linspl3)#(c.Meet_tavg c.Meet_tavg#c.Meet_tavg) $weather_control, absorb($absorb_var) vce(cluster meet_code athlete)
	estimates save "Semi_Parametric_Quadratic_Hetero", replace

*** linear spline with 2 knots 
	mkspline meet_linspl1 50 meet_linspl2 60 meet_linspl3 = Meet_tavg, di

	reghdfe norm_result c.meet_linspl* $weather_control, absorb($absorb_var) vce(cluster meet_code athlete)
	estimates save "Semi_Parametric_Linear_Spline_2Knots_Homo", replace

	reghdfe norm_result c.meet_linspl* (c.team_linspl1 c.team_linspl2 c.team_linspl3)#c.meet_linspl* $weather_control, absorb($absorb_var) vce(cluster meet_code athlete)
	estimates save "Semi_Parametric_Linear_Spline_2Knots_Hetero", replace

*** cubic spline with 3 knots
	mkspline2 meet_cubspl= Meet_tavg, cubic di knots(49.6 56.7 63.8)

	reghdfe norm_result c.meet_cubspl* $weather_control, absorb($absorb_var) vce(cluster meet_code athlete)
	estimates save "Semi_Parametric_Cubic_Spline_Homo", replace

	reghdfe norm_result c.meet_cubspl* (c.team_linspl1 c.team_linspl2 c.team_linspl3)#c.meet_cubspl* $weather_control, absorb($absorb_var) vce(cluster meet_code athlete)
	estimates save "Semi_Parametric_Cubic_Spline_Hetero", replace

*** cubic
	reghdfe norm_result Meet_tavg c.Meet_tavg#c.Meet_tavg c.Meet_tavg#c.Meet_tavg#c.Meet_tavg $weather_control, absorb($absorb_var) vce(cluster meet_code athlete)
	estimates save "Semi_Parametric_Cubic_Homo", replace

	reghdfe norm_result Meet_tavg c.Meet_tavg#c.Meet_tavg c.Meet_tavg#c.Meet_tavg#c.Meet_tavg (c.team_linspl1 c.team_linspl2 c.team_linspl3)#(c.Meet_tavg c.Meet_tavg#c.Meet_tavg c.Meet_tavg#c.Meet_tavg#c.Meet_tavg) $weather_control, absorb($absorb_var) vce(cluster meet_code athlete)
	estimates save "Semi_Parametric_Cubic_Hetero", replace

*-----------------------calculate performance changes--------------------------*
	
global RESULT "$PATH/Results/Simulations/Robustness/Min Max Avg"
cd "$PATH/Results/Simulations/Robustness/Min Max Avg"

foreach scenario in "RCP85" "RCP45" "1F" "2F" "4F" {

	*** load data set 
		use "$PATH/Data/TF_Data.dta", clear

	*** keep only endurance events
		keep if event_type == "Endurance"

	*** drop training temperature lower than 40F to mitigate noise
	/*avoid the case where athletes might get trained indoors during extreme cold weathers*/
		drop if Team_tavg <= 40 

	*** use average of tmin and tmax
		replace Meet_tavg = (Meet_tmax + Meet_tmin) / 2
		replace Team_tavg = (Team_tmax + Team_tmin_7) / 2

	*** set scenario
		if "`scenario'" == "RCP85" {
			gen Meet_GDDPtemp = Meet_RCP85temp
			gen Team_GDDPtemp = Team_RCP85temp 
		}
		if "`scenario'" == "RCP45" {
			gen Meet_GDDPtemp = Meet_RCP45temp
			gen Team_GDDPtemp = Team_RCP45temp 
		}
		if "`scenario'" == "1F" {
			gen Meet_GDDPtemp = Meet_tavg + 1
			gen Team_GDDPtemp = Team_tavg + 1
		}
		if "`scenario'" == "2F" {
			gen Meet_GDDPtemp = Meet_tavg + 2
			gen Team_GDDPtemp = Team_tavg + 2
		}
		if "`scenario'" == "4F" {
			gen Meet_GDDPtemp = Meet_tavg + 4
			gen Team_GDDPtemp = Team_tavg + 4
		}

	*** only keep weather variables
		keep athlete Meet_tavg Team_tavg Meet_GDDPtemp Team_GDDPtemp Team_zip norm_result

	*** set other control variables at the mean level
		gen wind = -0.0000445
		gen Meet_preciptotal = 0.0819
	    gen Meet_ozone = 0.0367131
	    gen travel_dist = 277.2341

	*** drop observations with missing data
		egen rowmiss = rowmiss(Meet_tavg Team_tavg Meet_GDDPtemp Team_GDDPtemp)
		drop if rowmiss > 0

	*** construct temperature variables
		* team temperature 
		mkspline team_c_linspl1 50 team_c_linspl2 60 team_c_linspl3 = Team_tavg, di
		mkspline team_f_linspl1 50 team_f_linspl2 60 team_f_linspl3 = Team_GDDPtemp, di

		* meet temperature - cubic spline
		mkspline2 meet_c_cubspl= Meet_tavg, cubic di knots(49.6 56.7 63.8)
		mkspline2 meet_f_cubspl= Meet_GDDPtemp, cubic di knots(49.6 56.7 63.8)

		* meet temperature - linear spline with 2 knots
		mkspline meet_c_linspl1 50 meet_c_linspl2 60 meet_c_linspl3 = Meet_tavg, di
		mkspline meet_f_linspl1 50 meet_f_linspl2 60 meet_f_linspl3 = Meet_GDDPtemp, di

		* meet temperature - linear spline with 3 knots
		mkspline meet_c_3linspl1 50 meet_c_3linspl2 60 meet_c_3linspl3 70 meet_c_3linspl4 = Meet_tavg, di
		mkspline meet_f_3linspl1 50 meet_f_3linspl2 60 meet_f_3linspl3 70 meet_f_3linspl4 = Meet_GDDPtemp, di

		* backup meet and team temperature
		gen Meet_tavg_bkp = Meet_tavg 
		gen Team_tavg_bkp = Team_tavg 

	*** initialize variables to store linear splines of team pre-exposure temperatures
		forvalues i = 1/3 {
			gen team_linspl`i' = .
		}

	*** linear spline with 3 knots
		forvalues i = 1/4 {
			gen meet_3linspl`i' = .
		}

		* scenario 1
		estimates use "Semi_Parametric_Linear_Spline_3knots_Homo"

		forvalues i = 1/4 {
			replace meet_3linspl`i' = meet_c_3linspl`i'
		}
		predict yhat_c1_linear_spline3, xb

		forvalues i = 1/4 {
			replace meet_3linspl`i' = meet_f_3linspl`i'
		}
		predict yhat_f1_linear_spline3, xb

		* scenario 2
		estimates use "Semi_Parametric_Linear_Spline_3knots_Hetero"

		forvalues i = 1/4 {
			replace meet_3linspl`i' = meet_c_3linspl`i'
		}
		forvalues i = 1/3 {
			replace team_linspl`i' = team_c_linspl`i'
		}
		predict yhat_c2_linear_spline3, xb

		forvalues i = 1/4 {
			replace meet_3linspl`i' = meet_f_3linspl`i'
		}
		forvalues i = 1/3 {
			replace team_linspl`i' = team_f_linspl`i'
		}
		predict yhat_f2_linear_spline3, xb


	*** quadratic
		* scenario 1
		estimates use "Semi_Parametric_Quadratic_Homo"

		replace Meet_tavg = Meet_tavg_bkp
		predict yhat_c1_quadratic, xb

		replace Meet_tavg = Meet_GDDPtemp
		predict yhat_f1_quadratic, xb

		* scenario 2
		estimates use "Semi_Parametric_Quadratic_Hetero"
		replace Meet_tavg = Meet_tavg_bkp
		forvalues i = 1/3 {
			replace team_linspl`i' = team_c_linspl`i'
		}
		predict yhat_c2_quadratic, xb

		replace Meet_tavg = Meet_GDDPtemp
		forvalues i = 1/3 {
			replace team_linspl`i' = team_f_linspl`i'
		}
		predict yhat_f2_quadratic, xb

	*** cubic
		* scenario 1
		estimates use "Semi_Parametric_Cubic_Homo"

		replace Meet_tavg = Meet_tavg_bkp
		predict yhat_c1_cubic, xb

		replace Meet_tavg = Meet_GDDPtemp
		predict yhat_f1_cubic, xb

		* scenario 2
		estimates use "Semi_Parametric_Cubic_Hetero"
		replace Meet_tavg = Meet_tavg_bkp
		forvalues i = 1/3 {
			replace team_linspl`i' = team_c_linspl`i'
		}
		predict yhat_c2_cubic, xb

		replace Meet_tavg = Meet_GDDPtemp
		forvalues i = 1/3 {
			replace team_linspl`i' = team_f_linspl`i'
		}
		predict yhat_f2_cubic, xb

	*** linear spline with 2 knots
		forvalues i = 1/3 {
			gen meet_linspl`i' = .
		}

		* scenario 1
		estimates use "Semi_Parametric_Linear_Spline_2Knots_Homo"
		forvalues i = 1/3 {
			replace meet_linspl`i' = meet_c_linspl`i'
		}
		predict yhat_c1_linear_spline2, xb

		forvalues i = 1/3 {
			replace meet_linspl`i' = meet_f_linspl`i'
		}
		predict yhat_f1_linear_spline2, xb

		* scenario 2
		estimates use "Semi_Parametric_Linear_Spline_2Knots_Hetero"
		forvalues i = 1/3 {
			replace meet_linspl`i' = meet_c_linspl`i'
		}
		forvalues i = 1/3 {
			replace team_linspl`i' = team_c_linspl`i'
		}
		predict yhat_c2_linear_spline2, xb

		forvalues i = 1/3 {
			replace meet_linspl`i' = meet_f_linspl`i'
		}
		forvalues i = 1/3 {
			replace team_linspl`i' = team_f_linspl`i'
		}
		predict yhat_f2_linear_spline2, xb

	*** cubic spline
		forvalues i = 1/2 {
			gen meet_cubspl`i' = .
		}

		* scenario 1
		estimates use "Semi_Parametric_Cubic_Spline_Homo"
		forvalues i = 1/2 {
			replace meet_cubspl`i' = meet_c_cubspl`i'
		}
		predict yhat_c1_cubic_spline, xb

		forvalues i = 1/2 {
			replace meet_cubspl`i' = meet_f_cubspl`i'
		}
		predict yhat_f1_cubic_spline, xb

		* scenario 2
		estimates use "Semi_Parametric_Cubic_Spline_Hetero"
		forvalues i = 1/2 {
			replace meet_cubspl`i' = meet_c_cubspl`i'
		}
		forvalues i = 1/3 {
			replace team_linspl`i' = team_c_linspl`i'
		}
		predict yhat_c2_cubic_spline, xb

		forvalues i = 1/2 {
			replace meet_cubspl`i' = meet_f_cubspl`i'
		}
		forvalues i = 1/3 {
			replace team_linspl`i' = team_f_linspl`i'
		}
		predict yhat_f2_cubic_spline, xb

	*** calculate percentage changes
		foreach model in linear_spline3 quadratic cubic linear_spline2 cubic_spline {
			gen percent1_`model' = (yhat_f1_`model' - yhat_c1_`model') / abs(norm_result)
			gen percent2_`model' = (yhat_f2_`model' - yhat_c2_`model') / abs(norm_result)
		}

	*** calculate mean changes
		collapse (mean) percent1* percent2*, by(ath Team_zip Team_tavg)
		foreach model in linear_spline3 quadratic cubic linear_spline2 cubic_spline {
		forvalues i = 1/2 {
			replace percent`i'_`model' = percent`i'_`model' * 100
		}
		}

		collapse (mean)	percent1_linear_spline3 percent1_cubic_spline percent1_linear_spline2 percent1_quadratic percent1_cubic ///
						percent2_linear_spline3 percent2_cubic_spline percent2_linear_spline2 percent2_quadratic percent2_cubic ///
				 (semean) se1percent_linear_spline3 = percent1_linear_spline3 ///
				 		  se1percent_cubic_spline   = percent1_cubic_spline ///
				 		  se1percent_linear_spline2 = percent1_linear_spline2 ///
				 		  se1percent_quadratic 	 	= percent1_quadratic ///
				 		  se1percent_cubic  		= percent1_cubic ///
				 		  se2percent_linear_spline3 = percent2_linear_spline3 ///
				 		  se2percent_cubic_spline   = percent2_cubic_spline ///
				 		  se2percent_linear_spline2 = percent2_linear_spline2 ///
				 		  se2percent_quadratic      = percent2_quadratic ///
				 		  se2percent_cubic          = percent2_cubic 
		gen i = _n
		reshape long percent1 percent2 se1percent se2percent, i(i) j(model) string
		order percent1 se1percent percent2 se2percent
		replace i = 2 if model == "_quadratic"
		replace i = 3 if model == "_cubic"
		replace i = 4 if model == "_linear_spline2"
		replace i = 5 if model == "_cubic_spline"
		sort i

	*** export results 
		export excel using "$RESULT/Estimated_Performance_Change_`scenario'.xlsx", firstrow(var) replace

}


********************************************************************************
* Table A8: estimated performance change using Interpolated daily average 
*    	    temperature projection
********************************************************************************

global RESULT "$PATH/Results/Simulations/Robustness/Daily Interpolation"
cd "$PATH/Results/Simulations/Baseline"

foreach scenario in "RCP85" "RCP45" {

	*** load data set 
		use "$PATH/Data/TF_Data.dta", clear

	*** keep only endurance events
		keep if event_type == "Endurance"

	*** drop training temperature lower than 40F to mitigate noise
		/*avoid the case where athletes might get trained indoors during extreme cold weathers*/
		drop if Team_tavg <= 40 

	*** set scenario
		if "`scenario'" == "RCP85" {
			gen Meet_GDDPtemp = Meet_RCP85tempavg
			gen Team_GDDPtemp = Team_RCP85tempavg 
		}
		if "`scenario'" == "RCP45" {
			gen Meet_GDDPtemp = Meet_RCP45tempavg
			gen Team_GDDPtemp = Team_RCP45tempavg 
		}
	
	*** only keep weather variables
		keep athlete Meet_tavg Team_tavg Meet_GDDPtemp Team_GDDPtemp Team_zip norm_result

	*** set other control variables at the mean level
		gen wind = -0.0000445
		gen Meet_preciptotal = 0.0819
	    gen Meet_ozone = 0.0367131
	    gen travel_dist = 277.2341

	*** drop observations with missing data
		egen rowmiss = rowmiss(Meet_tavg Team_tavg Meet_GDDPtemp Team_GDDPtemp)
		drop if rowmiss > 0

	*** construct temperature variables
		* team temperature 
		mkspline team_c_linspl1 50 team_c_linspl2 60 team_c_linspl3 = Team_tavg, di
		mkspline team_f_linspl1 50 team_f_linspl2 60 team_f_linspl3 = Team_GDDPtemp, di

		* meet temperature - cubic spline
		mkspline2 meet_c_cubspl= Meet_tavg, cubic di knots(49.6 56.7 63.8)
		mkspline2 meet_f_cubspl= Meet_GDDPtemp, cubic di knots(49.6 56.7 63.8)

		* meet temperature - linear spline with 2 knots
		mkspline meet_c_linspl1 50 meet_c_linspl2 60 meet_c_linspl3 = Meet_tavg, di
		mkspline meet_f_linspl1 50 meet_f_linspl2 60 meet_f_linspl3 = Meet_GDDPtemp, di

		* meet temperature - linear spline with 3 knots
		mkspline meet_c_3linspl1 50 meet_c_3linspl2 60 meet_c_3linspl3 70 meet_c_3linspl4 = Meet_tavg, di
		mkspline meet_f_3linspl1 50 meet_f_3linspl2 60 meet_f_3linspl3 70 meet_f_3linspl4 = Meet_GDDPtemp, di

		* backup meet and team temperature
		gen Meet_tavg_bkp = Meet_tavg 
		gen Team_tavg_bkp = Team_tavg 

	*** initialize variables to store linear splines of team pre-exposure temperatures
		forvalues i = 1/3 {
			gen team_linspl`i' = .
		}

	*** linear spline with 3 knots
		forvalues i = 1/4 {
			gen meet_3linspl`i' = .
		}

		* scenario 1
		estimates use "Semi_Parametric_Linear_Spline_3knots_Homo"

		forvalues i = 1/4 {
			replace meet_3linspl`i' = meet_c_3linspl`i'
		}
		predict yhat_c1_linear_spline3, xb

		forvalues i = 1/4 {
			replace meet_3linspl`i' = meet_f_3linspl`i'
		}
		predict yhat_f1_linear_spline3, xb

		* scenario 2
		estimates use "Semi_Parametric_Linear_Spline_3knots_Hetero"

		forvalues i = 1/4 {
			replace meet_3linspl`i' = meet_c_3linspl`i'
		}
		forvalues i = 1/3 {
			replace team_linspl`i' = team_c_linspl`i'
		}
		predict yhat_c2_linear_spline3, xb

		forvalues i = 1/4 {
			replace meet_3linspl`i' = meet_f_3linspl`i'
		}
		forvalues i = 1/3 {
			replace team_linspl`i' = team_f_linspl`i'
		}
		predict yhat_f2_linear_spline3, xb


	*** quadratic
		* scenario 1
		estimates use "Semi_Parametric_Quadratic_Homo"

		replace Meet_tavg = Meet_tavg_bkp
		predict yhat_c1_quadratic, xb

		replace Meet_tavg = Meet_GDDPtemp
		predict yhat_f1_quadratic, xb

		* scenario 2
		estimates use "Semi_Parametric_Quadratic_Hetero"
		replace Meet_tavg = Meet_tavg_bkp
		forvalues i = 1/3 {
			replace team_linspl`i' = team_c_linspl`i'
		}
		predict yhat_c2_quadratic, xb

		replace Meet_tavg = Meet_GDDPtemp
		forvalues i = 1/3 {
			replace team_linspl`i' = team_f_linspl`i'
		}
		predict yhat_f2_quadratic, xb

	*** cubic
		* scenario 1
		estimates use "Semi_Parametric_Cubic_Homo"

		replace Meet_tavg = Meet_tavg_bkp
		predict yhat_c1_cubic, xb

		replace Meet_tavg = Meet_GDDPtemp
		predict yhat_f1_cubic, xb

		* scenario 2
		estimates use "Semi_Parametric_Cubic_Hetero"
		replace Meet_tavg = Meet_tavg_bkp
		forvalues i = 1/3 {
			replace team_linspl`i' = team_c_linspl`i'
		}
		predict yhat_c2_cubic, xb

		replace Meet_tavg = Meet_GDDPtemp
		forvalues i = 1/3 {
			replace team_linspl`i' = team_f_linspl`i'
		}
		predict yhat_f2_cubic, xb

	*** linear spline with 2 knots
		forvalues i = 1/3 {
			gen meet_linspl`i' = .
		}

		* scenario 1
		estimates use "Semi_Parametric_Linear_Spline_2Knots_Homo"
		forvalues i = 1/3 {
			replace meet_linspl`i' = meet_c_linspl`i'
		}
		predict yhat_c1_linear_spline2, xb

		forvalues i = 1/3 {
			replace meet_linspl`i' = meet_f_linspl`i'
		}
		predict yhat_f1_linear_spline2, xb

		* scenario 2
		estimates use "Semi_Parametric_Linear_Spline_2Knots_Hetero"
		forvalues i = 1/3 {
			replace meet_linspl`i' = meet_c_linspl`i'
		}
		forvalues i = 1/3 {
			replace team_linspl`i' = team_c_linspl`i'
		}
		predict yhat_c2_linear_spline2, xb

		forvalues i = 1/3 {
			replace meet_linspl`i' = meet_f_linspl`i'
		}
		forvalues i = 1/3 {
			replace team_linspl`i' = team_f_linspl`i'
		}
		predict yhat_f2_linear_spline2, xb

	*** cubic spline
		forvalues i = 1/2 {
			gen meet_cubspl`i' = .
		}

		* scenario 1
		estimates use "Semi_Parametric_Cubic_Spline_Homo"
		forvalues i = 1/2 {
			replace meet_cubspl`i' = meet_c_cubspl`i'
		}
		predict yhat_c1_cubic_spline, xb

		forvalues i = 1/2 {
			replace meet_cubspl`i' = meet_f_cubspl`i'
		}
		predict yhat_f1_cubic_spline, xb

		* scenario 2
		estimates use "Semi_Parametric_Cubic_Spline_Hetero"
		forvalues i = 1/2 {
			replace meet_cubspl`i' = meet_c_cubspl`i'
		}
		forvalues i = 1/3 {
			replace team_linspl`i' = team_c_linspl`i'
		}
		predict yhat_c2_cubic_spline, xb

		forvalues i = 1/2 {
			replace meet_cubspl`i' = meet_f_cubspl`i'
		}
		forvalues i = 1/3 {
			replace team_linspl`i' = team_f_linspl`i'
		}
		predict yhat_f2_cubic_spline, xb

	*** calculate percentage changes
		foreach model in linear_spline3 quadratic cubic linear_spline2 cubic_spline {
			gen percent1_`model' = (yhat_f1_`model' - yhat_c1_`model') / abs(norm_result)
			gen percent2_`model' = (yhat_f2_`model' - yhat_c2_`model') / abs(norm_result)
		}

	*** calculate mean changes
		collapse (mean) percent1* percent2*, by(ath Team_zip Team_tavg)
		foreach model in linear_spline3 quadratic cubic linear_spline2 cubic_spline {
		forvalues i = 1/2 {
			replace percent`i'_`model' = percent`i'_`model' * 100
		}
		}

		collapse (mean)	percent1_linear_spline3 percent1_cubic_spline percent1_linear_spline2 percent1_quadratic percent1_cubic ///
						percent2_linear_spline3 percent2_cubic_spline percent2_linear_spline2 percent2_quadratic percent2_cubic ///
				 (semean) se1percent_linear_spline3 = percent1_linear_spline3 ///
				 		  se1percent_cubic_spline   = percent1_cubic_spline ///
				 		  se1percent_linear_spline2 = percent1_linear_spline2 ///
				 		  se1percent_quadratic 	 	= percent1_quadratic ///
				 		  se1percent_cubic  		= percent1_cubic ///
				 		  se2percent_linear_spline3 = percent2_linear_spline3 ///
				 		  se2percent_cubic_spline   = percent2_cubic_spline ///
				 		  se2percent_linear_spline2 = percent2_linear_spline2 ///
				 		  se2percent_quadratic      = percent2_quadratic ///
				 		  se2percent_cubic          = percent2_cubic 
		gen i = _n
		reshape long percent1 percent2 se1percent se2percent, i(i) j(model) string
		order percent1 se1percent percent2 se2percent
		replace i = 2 if model == "_quadratic"
		replace i = 3 if model == "_cubic"
		replace i = 4 if model == "_linear_spline2"
		replace i = 5 if model == "_cubic_spline"
		sort i

	*** export results 
		export excel using "$RESULT/Estimated_Performance_Change_`scenario'.xlsx", firstrow(var) replace

}


********************************************************************************
* Figure A1: distribution of competition and training temperatures
******************************************************************************** 

*** set results folder
	cd "$PATH/Results/Simulations/Baseline"

*** load data set 
	use "$PATH/Data/TF_Data.dta", clear

*** keep only endurance events
	keep if event_type == "Endurance"

*** keep variables of interest
	keep Meet_tavg Team_tavg Meet_RCP* Team_RCP* athlete Team_zip Meet_zip 

*** plotting
	grstyle init 
	grstyle set plain, nogrid compact

	* competition temperatures
	histogram Meet_tavg ///
		, bin(45) fcolor(gs12) lcolor(gs14) xlabel(35(5)85) xtitle(Temperature (F)) fxsize(50) ///
		title(Historical) ylabel(0(0.02)0.06) ///
		saving("Distribution_Meet_Tavg.gph", replace)
	histogram Meet_RCP45temp ///
		, bin(45) fcolor(gs12) lcolor(gs14) xlabel(35(5)85) xtitle(Temperature (F)) fxsize(50) ///
		title(RCP4.5) ylabel(0(0.02)0.06) ///
		saving("Distribution_Meet_RCP45.gph", replace)
	histogram Meet_RCP85temp ///
		, bin(45) fcolor(gs12) lcolor(gs14) xlabel(35(5)85) xtitle(Temperature (F)) fxsize(50) ///
		title(RCP8.5) ylabel(0(0.02)0.06) ///
		saving("Distribution_Meet_RCP85.gph", replace)
	graph combine ///
		"Distribution_Meet_Tavg.gph" "Distribution_Meet_RCP45.gph" "Distribution_Meet_RCP85.gph" ///
		, col(1) imargin(0 0.6 0 0) xsize(10) ysize(12) scale(1) xcommon title("     Competition")
	graph save "Distribution_Meet.gph", replace
		
	* training temperatures
	histogram Team_tavg ///
		, bin(45) fcolor(gs12) lcolor(gs14) xlabel(35(5)85) xtitle(Temperature (F)) fxsize(50) ///
		title(Historical) ylabel(0(0.02)0.06) ///
		saving("Distribution_Team_Tavg.gph", replace)
	histogram Team_RCP45temp ///
		, bin(45) fcolor(gs12) lcolor(gs14) xlabel(35(5)85) xtitle(Temperature (F)) fxsize(50) ///
		title(RCP4.5) ylabel(0(0.02)0.06) ///
		saving("Distribution_Team_RCP45.gph", replace)
	histogram Team_RCP85temp ///
		, bin(45) fcolor(gs12) lcolor(gs14) xlabel(35(5)85) xtitle(Temperature (F)) fxsize(50) ///
		title(RCP8.5) ylabel(0(0.02)0.06) ///
		saving("Distribution_Team_RCP85.gph", replace)
	graph combine ///
		"Distribution_Team_Tavg.gph" "Distribution_Team_RCP45.gph" "Distribution_Team_RCP85.gph" ///
		, col(1) imargin(0 0.6 0 0) xsize(10) ysize(12) scale(1) xcommon title("     Training")
	graph save "Distribution_Team.gph", replace

	* combine figures
	graph combine ///
		"Distribution_Meet.gph" "Distribution_Team.gph" ///
		, col(2) imargin(0 0.6 0 0) xsize(10) ysize(12) scale(1) xcommon
	graph export "Temperature_Distribution.eps", replace

*** delete intermediate figures 
	foreach v in Meet Team {
		foreach x in Tavg RCP45 RCP85 {
			erase "Distribution_`v'_`x'.gph"
		}
		erase "Distribution_`v'.gph"
	}
